#####
#####
###     Estimate nonmetric vector model of preferences 
###     for seven values, using ranked value choices
###     from post-election wave of 2006 CCES, plot
###     estimated vector model, calculate and plot
###     mean vector, and write out text file containing
###     CCES case ID and vector coordinates for
###     each respondent.
###
###     Originally run August 16, 2011
###
###     Archived version prepared June 23, 2014
###
###     Note that this R session requires the 
###     "optiscale" package to carry out the
###     optimal scaling.
#####
#####

library(optiscale)
library(lattice)

###
###   Read data from file, "ranked
###   values, with names.txt"
###

vals <- read.table(file.choose(), header = TRUE)

values <- vals[,2:8]

values[1:10,]

###
###   Standardize data within rows
###

vals.std <- t(apply(values, 1, scale))

vals.std[1:10,]

###
###   Initialize the matrix of optimally-scaled
###   values, the previous fit, the iteration
###   number, and the improvement in fit over
###   the previous iteration
###

vals.os <- vals.std

prev.fit2 <- 0

niter <- 0

improve = 1

###
###   Start iterations
###

   while (improve > .01 & niter <= 25) {
   niter <- niter + 1

###
###   Perform SVD on OS version of rank-orders
###

   decomp <- svd(vals.os)

###
###   Calculate 2-dimensional goodness-of-fit
###   and improvement in fit over previous iteration
###
   
   d.sqd <- decomp$d^2
   
   fit2 <- sum(d.sqd[1:2]) / sum(d.sqd)
   
   fitvector <- cumsum(d.sqd) / rep(sum(d.sqd), length(d.sqd))
   
   improve <- fit2 - prev.fit2

###
###   Create iteration history
###
   
   if (niter == 1) {
      history <- c(niter, fit2, improve)
      }
   if (niter > 1)  {
      history <- rbind(history, c(niter, fit2, improve))
      }

###
###   Obtain terminal points of vectors for row objects
###   (respondents in this case), calculate predicted
###   ranks, norm the row object vectors to unit length
###
  
   subj.c1 <- decomp$u[, 1:2] %*% diag(decomp$d[1:2])
   
   pred.vals <- subj.c1 %*% t(decomp$v[,1:2])

   sumsqd <- apply(subj.c1^2, 1, sum)

   root.sumsqd <- (matrix(sumsqd, nrow = length(sumsqd), ncol = 1)) ^ .5
   
   subj.coord <- subj.c1 / (root.sumsqd %*% matrix(1, nrow = 1, ncol = 2))

###
###   Obtain new optimally scaled data values,
###   using predicted ranks from fitted model
###

   for (i in 1:nrow(vals.os)) {
      opped <- opscale(x.qual = vals.std[i,],
                             x.quant = pred.vals[i,],
                             level = 2,
                             process = 1,
                             rescale = T)
      vals.os[i,] <- opped$os
      }

###
###   Set current fit to previous fit for next iteration
###
   
   prev.fit2 <- fit2
   
   }

###
###   Print iteration history
###

history

###
###   The "fitvector" object shows the R-squared
###   for the solution in each dimensionality
###

fitvector

###
###   Plot value points
###

val.coords <- as.data.frame(decomp$v[, 1:2])

rownames(val.coords) <- colnames(values)

val.coords

xyplot(V2 ~ V1, val.coords,
   aspect = 1,
   xlim = c(-.8, .8),
   ylim = c(-.8, .8),
   panel = function (x, y) {
   panel.xyplot(x, y, col = "black")
   panel.text(x, y, labels = rownames(val.coords),
      pos = 1, cex = .75)
   }
)

###
###   Plot terminal points of row object vectors
###   in same space as value points
###

subj.coord <- as.data.frame(subj.coord)

set.seed(123)
xyplot(V2 ~ V1, val.coords,
   aspect = 1,
   xlim = c(-1.2, 1.2),
   ylim = c(-1.2, 1.2),
   panel = function (x, y) {
   panel.xyplot(x, y, col = "black", pch = 16)
   panel.xyplot(jitter(subj.coord$V1, amount = .05), 
      jitter(subj.coord$V2, amount = .05), col = "black")
   panel.text(x, y, labels = rownames(val.coords),
      pos = 1, cex = .75)
   }
)

###
###   Calculate overall mean direction and
###   mean resultant length
###

meand1 <- mean(subj.coord$V1)

meand2 <- mean(subj.coord$V2)

meand1 

meand2

mean.length <- (meand1^2 + meand2^2)^.5

mean.length

###
###   Plot mean vector along with 
###   individual vectors and value points
###

set.seed(123)
xyplot(V2 ~ V1, val.coords,
   aspect = 1,
   xlim = c(-1.2, 1.2),
   ylim = c(-1.2, 1.2),
   panel = function (x, y) {
   panel.xyplot(x, y, col = "black", pch = 16)
   panel.xyplot(jitter(subj.coord$V1, amount = .05), 
      jitter(subj.coord$V2, amount = .05), col = "black")
   panel.text(x, y, labels = rownames(val.coords),
      pos = 1, cex = .75)
   panel.arrows(0, 0, meand1, meand2, angle = 15, length = .15)
   }
)


###
###   Merge data on row object vectors with case IDs from
###   original dataset
###

outvecs <- cbind(vals[,1], subj.coord)

outvecs[, 2:3] <- round(outvecs[, 2:3], digits = 5)

###
###   Write data on row object vectors to external file,
###   for merging back in to CCES data
###

write.table(outvecs, file = "row vectors.txt",
   sep = " ", 
      row.names = FALSE,
      col.names = FALSE
)
